Nonlinear optics of III-V semiconductors in the terahertz regime: an ab-initio study 
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We compute from first principles the infrared dispersion of the nonlinear susceptibility x in 
zincblende semiconductors. At terahertz frequencies the nonlinear susceptibility depends not only 
on the purely electronic response Xc»' , but also on three other parameters Ci , C2 and C3 describing 
the contributions from ionic motion. They relate to the TO Raman polarizability, the second- 
order displacement-induced dielectric polarization, and the third-order lattice potential. Contrary to 
previous theory, we find that mechanical anharmonicity (C3) dominates over electrical anharmonicity 
(C2), which is consistent with recent experiments on GaAs. We predict that the sharp minimum in 
the intensity of second-harmonic generation recently observed for GaAs between wto/S and wto 
does not occur for several other III-V compounds. 

PACS numbers: 78.30.Fs, 71.36.+C, 78.20.Jq, 42.65.An 
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I. INTRODUCTION 

The nonhnear optical properties of materials in the vis- 
ible and near-infrared (IR) have been extensively studied 
since the development of the laser. Except for a few pi- 
oneering effortsjiiSi the far-IR region of the spectrum has 
remained largely unexplored. This state of affairs, a con- 
sequence of the lack of tunable and intense laser sources 
and sensitive detectors in the terahertz range, is start- 
ing to change, thanks to advances in instrumentation,?'^ 
More accurate measurements of the nonlinear suscepti- 
bilities at terahertz frequencies are beginning to appear,? 
calling for quantitative theoretical modeling. 



The first nonlinear susceptibility x''^'' (^3 



(^2; wi, a;2), which is nonzero in acentric crystals, displays 
strong dispersion when the frequencies involved are near 
a zone-center transverse optical phonon frequency wto- 
This behavior was first observed by Faust and Henry^ 
(see also Refs.j^andQ) in GaP for the case of mixing be- 
tween visible and far-IR radiation (wi, o^i + L02 > t^TO > 
UI2)] they showed that the dispersion of this process de- 

pends on the linear electro-optic susceptibility Xeo • An- 
other early set of experimentsSiS investigated frequency 
mixing in the microwave range below the lattice reso- 
nances (wi,ti;2,wi + UJ2 < i^To)- For the two zincblende 
compounds studied, GaAs and GaP, it was found^ that 
the sign of the microwave coefficient Xmw (the static non- 
linear susceptibility) was opposite to that of the "high- 

frequency" coefficient Xoo describing second-harmonic 
generation (SHG) and frequency mixing in the trans- 
parency region of the crystal (wto < ^1,^2,^1 + ^2 < 
Eg/Ti). Mayer and Keilmann^ later studied the disper- 
sion of the SHG coefficient Xshg('^) ~ X^'^\'^^]^^^) of 
GaAs and LiTaOs over a limited frequency range (0.6 to 
1.7 THz). 

The present work was motivated by the recent experi- 
ments of Dekorsy et al. on SHG in GaAsn^- Using a tun- 



able free-electron laser they measured the dispersion of 

(2) 

XsHG from 4 to 6 THz, observing a strong resonant en- 
hancement at 4.5 THz, close to lo'yo/2 ~ 4.1 THz as 
expected, followed by a sharp dip in the power output. 

Although they were unable to determine unambiguously 

(2) 

the frequency cjo at which Xshg vanishes, a lower bound 
of 5.3 THz was established. By fitting the location of 
the minimum in SHG power to an expression derived by 
Flytzanis,^° they obtained new parameters for GaAs. 

Flytzanis's expression is given below [Eq. (Q]. Us- 
ing an effective-bond model, he obtained numerical es- 
timates for its parameters. These have been used as a 
guide for interpreting subsequent experiments?^ in spite 
of their questionable reliability. For instance, the model 
predicts the wrong sign for the Born effective charge and 
the Raman polarizability, two of the material parameters 
at play. In this work we re-examine this problem using 
first-principles density-functional techniques. 

The paper is organized as follows. In Sec. |n] we dis- 
cuss the formalism describing the IR dispersion of x'^''- 
The computational approach is explained in Sec. Illll In 
Sec. IIVI we present and discuss results for several III-V 
semiconductors, GaAs, GaP, AlP, AlAs, and AlSb. Fi- 
nally, in Sec. we summarize our findings. Additional 
discussion of the computational methods is given in two 
Appendices. 



II. FORMALISM 

In this work we limit ourselves to the zincblende struc- 
ture adopted by III-V semiconductors, the simplest crys- 
tal structure where a non-vanishing x^^^ is allowed by 
symmetry. In zincblende there is a single TO mode, and 
third-rank tensors such as x''^'' have only one independent 
component x^xyz- The dispersion of x'^'' below the elec- 
tronic resonances and above the elastic resonances of the 
medium is given by the following expression, obtained by 
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Flytzanis>i£ 

X^^H^i +W2;wi,w2) = A(^i,'^2), (1) 



where 

k(uJi,UJ2) = 1 + Ci 

+C3 



1 



1 



1 



^(Wi) D{U02) D{lJi+LU2) 
1 1 



i:i(wi)D(w2) D{UJ2)D{UJI+LU2) 
1 



1 

D{LJi)D{Lj2)D{iL;i +UJ2) 



(2) 



and D(uj) = 1 — cj^/cj^q — ijuj/to^Q is the resonance de- 
nominator with phonon damping. This expression can 
be derived under rather general conditions^ and is in- 
dependent of the details of the microscopic forces, the 
only assumption being that the mode is weakly damped 

(7 < UJTo)- 

The dispersion depends on three dimensionless coeffi- 
cients, 



aTO ( 



2wx 



(2) 



(known as the Faust-Henry coefficient), 



C2 



,(2) 



2vx 



(2) 



and 



C. = - 



0(3) 



2vx 



(2) 



z* 



(3) 



(4) 



(5) 



Here v is the volume of the primitive cell, M is the re- 
duced mass, and the remaining quantities are discussed 
below. Unlike the second-rank tensor properties, the 
signs of the C^'s remain unchanged if we reverse the 
definition of the positive [111] direction. We adopt the 

convention that it points from the cation to the closest 
anioni^2,i3 

Having obtained the above model-independent expres- 
sion for the IR dispersion of x*"^"*) Flytzanis then pro- 
ceeded to estimate the values of the coefficients Ci, C2, 
and C3 for several HI-V compounds, using an effective- 
bond model. We will now describe how to compute 
them from first-principles. Except for the damping pa- 
rameter 7 (which we do not calculate), the quantities 
entering Eqs. are conveniently evaluated as deriva- 

tives of forces or macroscopic polarization with respect 
to macroscopic electric fields or displacements (forces un- 
der finite fields are readily available via the Hellmann- 
Feynman theorem—). In what follows F represents the 
force on the cation, P is the macroscopic polarization, S 
is a macroscopic electric field, and u = uni — uy is the 



relative displacement between the cation (group III) and 
anion (group V) sublattices away from their equilibrium 
positions. Two of the quantities in Eqs. (mSJ are first 
derivatives: the cation Born effective charge 



Z* 



E=0 



d£x 



u=0 



and the zone-center TO phonon frequency 



M dur 



(6) 



(7) 



£=0 



The remaining four are second derivatives: the non- 
resonant electronic ("high-frequency") quadratic suscep- 
tibility 



.(2) _ 



1 d^P. 



2 dEM,, 



(8) 



u=0 



the non-resonant TO Raman polarizability per primitive 
cell^ 



QfTO 



d[x. 



(1)1 

00 \xy 



di 



£=0 



d£xd£y 



(9) 



u=0 



the second-order dipole moment or "electrical anhar- 
monicity" 



,(2) 



duzdu^j 



(10) 



£=0 



and the third-order lattice potential or "mechanical an- 
harmonicity" 



a(3) - 



d^Fx 



duzduy 



(11) 



£=0 



Note that in Eq. ^ the Raman tensor was recast as a 
second-order field-induced forceii^^i^ 

It will be useful to consider Eq. JQ) in the three lim- 
iting cases discussed in the Introduction. In the high- 
frequency limit it reduces to the purely electronic coef- 
(2) 

ficient Xoo ; for uji,uji + uj2 > ujto > ^2 it becomes the 
unclamped-ion, strain-free electro-optic coefficient, 

=X^^\l + Ci); (12) 

and finally for wi, a;2, ti-"! + 1^2 < 'j-'to it describes the 
strain-free static (microwave) nonlinear susceptibility, 
which involves all three coefficients. 



^(2) 



X^^\l + 3C,+3C2 + C3). 



(13) 



III. COMPUTATIONAL APPROACH 



The calculations are done using ABINIT^ a plane- 
wave pseudopotential density-functional code, using both 
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the local-density approximation (LDA) and a general- 
ized gradient approximatiorJ^ (GGA). In order to bet- 
ter assess the sensitivity of the nonlinear optical prop- 
erties to the approximate density functional used, all 
calculations are performed at the same (experimental) 
lattice constants. Norm-conserving TrouUier-Martins 
pseudopotentials^^ are used for all materials, and for Ga 
the d-electrons are included in the valence. We used a 
cutoff energy of 30 Ha for the aluminum compounds, and 
45 Ha for the gallium compounds. 

For finite systems such as molecules, ab-initio calcula- 
tions of x^^) below the electronic resonances (including 
contributions from ionic motion2£) are performed rou- 
tinely. Similar calculations for bulk solids have become 
feasible only recently, thanks to a series of developments, 
starting with the Berry-phase theory of polarization,^ 
which provides a means of evaluating the electronic con- 
tribution to the macroscopic polarization P as a bulk 
quantity (for a review, see Ref. l22f) . In Ref. a density- 
functional perturbation method was developed for com- 
puting nonlinear (electronic, electro-optic, and nonreso- 
nant Raman) susceptibilities using the Berry-phase for- 
malism to treat the electric field perturbation. Here we 
use a closely-related approach where the derivatives are 
evaluated by finite differences rather than analytically, 
using the method of Ref. 14 to handle finite electric fields. 
The appeal of this method lies in its simplicity: once 
implemented, no extra coding is required to compute a 
different higher-order or mixed derivative, or to switch 
from, e.g., LDA to GGA. 

In our calculations we apply the finite perturbation (u 
or £) along [111], and monitor the response (P or F) 
along the same direction. That is, we let u = Sa{x + 
y + z), where a is the lattice constant and the cation 
and anion sublattices are brought closer together when 
d > 0] £ — e£'o(x + y + z), where Eq — e/(47reoaQ); 
F = f {± + y + z); and AP = p(x + y + z). Eqs. P JTT|) 
then become 



Z* 



a dp 
TdS 



dl 

de ' 



2 

^TO 



_J_dl 

aM 85 ' 



axo 



4 9e2' 



15V 
2 9e2' 



,(2) ^ ^ 



a d^p 



and 



8(9(52' 



2a2 86^ ■ 



(14) 
(15) 
(16) 
(17) 
(18) 

(19) 



The parameters are calculated from these expressions, 
and then inserted into Eqs. (PHHIl to obtain the coef- 
ficients Ci , C2 , and C3 . In Appendix El we describe 

a different approach whereby Xoo"* and Xmw are eval- 

(2) 

uated directly by finite differences in addition to Xoo , 
and Eqs. (|12|l and (|13|1 are then used to obtain Ci and 
3C2 + C3. 

We have taken as the smallest increments S = Ix 10~^ 
and e = 1 X IQ-* for AlAs, AlP, GaP, and AlSb. In the 
case of GaAs a smaller field step of £ = 3 x 10~^ was 
used. This was needed in order to stay below the A:-mesh 
dependent critical field above which the electric enthalpy 
functional loses its minima. Because of its smaller band 
gap, the critical field is lower for GaAs than for the other 
compounds. 

For the derivatives we use Richardson's extrapolation 
to estimate the limit h —> from calculations with two 
different step sizes: 

/(")(x) = ^i?("Hx,/i) - iz?(")(x,2/.) + 0(/i4), (20) 

where Z?^^-* is given by the centered finite difference ex- 
pression 

D^'^ (., h) . n-+'^)-n--'^) = ;^(,) + oih% 

(21) 

and Z?^2) is given by 



^(2)(^^/^) ^ fi^ + h) + fix ~h)^ 2f{x) 
= r{h) + 0{h^). 



(22) 



In order to speed up the convergence of polarization- 
dependent quantities with respect to the fc-point sam- 
pling, we use a similar extrapolation for the discretized 
Berry-phase formula, as described in Appendix All 
the values quoted in the tables were calculated on a 
16 X 16 X 16 fc-point mesh, except in the case of GaAs 
which, for reasons discussed in that Appendix, demanded 
a denser mesh. 



IV. RESULTS FOR III-V SEMICONDUCTORS 
A. Microscopic parameters 

We have systematically computed the values of all 
six parameters (j^ lll|l . using the methods summarized 
above, for five III-V compounds. The results are col- 
lected in Table Ifl and we will begin by discussing the 
comparison with the model theory of Flytzanis, and then 
pass to the discussion of the experimental values. 

Table |2 shows striking differences between the first- 
principles results and Flytzanis's model calculations: (i) 
While both levels of theory produce the same sign for 



(2) 



Xoo- 



they disagree on the signs of Z*, axO: /^^^\ and 
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TABLE I; Parameters that determine the nonhnear sus- 
ceptibihty of zinclende compounds in the infrared range, 
grouped into first- and second-derivative quantities [Eqs. JS]- 
lim l. Rows labelled "Model" pertain to the empirical- model 



of Flytzanis 



10.32 



calculations of Flytzanis: x 
maining values from Ref. 10 



(2) 



is taken from Ref. 124. the re- 







First derivative 


Second derivative 






Z* 




(2) 














(THz) 


(pm/V) (A^) 


(nC/m) 


(TJ/m^ 


GaAs LDA 


2.05 


8.0 


A 70 


— 04 


-1.89 


4.29 




GGA 


2.05 


8.3 


OO ( 


— 04t 


-1.67 


5.10 




Model 






127 


84 


5.90 


-1.2 




Expt. 


2.2" 


8.0" 


170* 








GaP 


LDA 


2.10 


10.6 


131 


-12.9 


-1.50 


4.82 




GGA 


2.15 


10.9 


114 


-9.8 


-1.36 


4.76 




Model 


-2.0 


11.5 


93 


43 


7.77 


-2.1 




Expt. 


2.0" 


11.0" 


71' 


-20" 






AlP 


LDA 


2.24 


12.7 


45 


-5 


0.28 


4.64 




GGA 


2.24 


13.1 


42 


-5 


0.24 


4.32 




Expt. 


2.28" 


13. 2"* 










AlAs 


LDA 


2.14 


10.4 


79 


-9 


0.27 


4.00 




GGA 


2.11 


10.8 


73 


-8 


0.22 


4.94 




Expt. 


2.3" 


10.8" 










AlSb 


LDA 


1.86 


9.2 


205 


-19 


-0.09 


3.09 




GGA 


1.79 


9.5 


187 


-18 


-0.12 


3.34 




Model 


-1.6 


9.8 


47 


77 


8.34 


-1.5 




Expt. 


1.9" 


9.6" 


153' 









"Quoted in Ref. .25, 
'Ref. HI 
"Ref.|23 

■^Quoted in Ref. ISl 
"Quoted in Ref. |2§ 



(ii) We find that the magnitude of fi^^^ ('Z''-^-') is 
significantly smaller (larger) than Flytzanis predicts. 

Born charges and phonon frequencies are routinely 
computed from first-principles, and they tend to com- 
pare favorably with experiment, as evidenced in Ta- 
ble Dielectric susceptibilities and Raman polarizabili- 
ties are more problematic. For example, it is well-known 
that density-functional theory tends to overestimate the 
dielectric constant. This also seems to be generally the 
case for the nonlinear susceptibility x'S ^Si^ The prob- 
lem here is compounded by the fact that the experimen- 
tal determination of this quantity is also problematic, 
with the values reported in the literature displaying a 
very large dispersionn^i We have opted for quoting the 
recommended values from Landolt-B6rnsteiniS& Inspec- 
tion of Table suggests that our first-principles values 
are too large by roughly a factor of two, which however 
is comparable with the uncertainty in the experimental 
determination. Measurements of the absolute Raman 
polarizability are also difficult, and few values are re- 
ported in the literature. Our result that axo < means 
that the bond polarizability along [111] increases with 
increasing bond length around the equilibrium length of 
the bond. This is in agreement with the measured sign in 
GaAs,^'^ but in disagreement with the model calculations 



(2) 



B. Lattice-induced contributions to x 

^From the calculated parameters in Table HI we ob- 
tained, using Eqs. (PHSJ, the various lattice-induced con- 
tributions to x^^\ collected in Table IITI We find that the 
Faust-Henry coefficients Ci are roughly a factor of two 
smaller than the experimental values, consistent with the 
overstimation of Xm"* in Eq. ^ discussed above. The me- 
chanical anharmonicity coefficient C3 is, in all cases, sig- 
nificantly larger in magnitude than C2 (electrical anhar- 
monicity). This is the opposite conclusion from Refs. ITol 
and We remark that although the correct signs for 
Ci and C3 were obtained therein, this resulted from a 
cancellation of errors in the signs of Z*, uto, and 0^^^ 
in Eqs. ^ and lO). No such cancellation occurs in C2, 
and indeed for the three compounds studied both in the 
present work and in Ref.^3 (AlSb, GaP and GaAs), there 
is a disagreement in the predicted sign for this quantity. 

While the coeffficicnt Ci can be measured in vari- 
ous ways (electro-optic effect, frequency mixing^ and 
relative Raman scattering efficiencies from LO and TO 
phonons'^'*'^''^) with fairly consistent results, it is difficult 
to disentangle the values of C2 and C3 from experiments. 
The more readily accessible quantity is the combination 
3C2 + C3: it follows from Eqs. ^ and ^ that 



3C2 + C3 = 2 



(2) 



(2) 
Xoo' 



(23) 



where all the quantities on the right-hand-side are di- 
rectly measurable. Flytzanis found C2 > and C3 < 
for all the III-V compounds he investigated. Under those 
circumstances, the measured sign of 3C2 + C3 indicates 
which anharmonic contribution (electrical or mechanical) 
is dominant in a given material. We find, however, that 
the sign of C2 is not the same for all III-V compounds, 
which invalidates such an analysis. 

Interestingly, our calculated 3C2 + C3 disagree in sign 
with the values inferred from experimenli^ (see Ta- 
ble ^Oj. It is apparent from Eq. (|23|l that the sign of 
3C2 + C3 is rather sensitive to not only the signs, but 

(2) (2) 

also the relative magnitudes of Xeo and Xmw- While the 
signs of our calculated Xmw and Xco"* agree with experi- 
ment (see Table InUl there is a significant discrepancy re- 
garding their magnitudes. The well-known limitations of 
density-functional theory in reproducing dielectric prop- 
erties, such as the optical gap (underestimated) and the 
the dielectric constant (overestimated) may be of concern 

in this regard. We note however that a possible error in 

(2) 

the magnitude of Xoo will not affect the sign of 3C2 + C3, 
while Z*, and the remaning parameters enter- 
ing C2 and C3, are expected to be reasonably accurate 
within density functional theory, which typically describe 
rather well lattice-dynamical and zero-field polarization 
properties (e.g.. Born charges) of III-V semiconductors!^ 
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TABLE II: Lattice- induced contributions to the nonlinear sus- 
ceptibility of zincblende compounds [Eqs. J^HSJ]. Rows la- 
belled "Model" pertain to the empirical-model calculations of 
Flytzanis. 







Ci 


C2 


C3 


AlP 


LDA 


-0.38 


0.05 


-1.82 




GGA 


-0.37 


0.04 


-1.78 


AlAs 


LDA 


-0.37 


0.03 


-0.91 




GGA 


-0.34 


0.02 


-0.82 


AlSb 


LDA 


-0.25 


-0.00 


-0.22 




GGA 


-0.23 


-0.00 


-0.18 




Model 


-1.97 


0.35 


-0.11 


GaP 


LDA 


-0.27 


-0.04 


-0.53 




GGA 


-0.28 


-0.07 


-0.56 




Model 


-0.37 


0.11 


-0.05 




Expt. 


—0.53" 








Expt. 


-0.75" 






GaAs 


LDA 


-0.35 


-0.02 


-0.12 




GGA 


-0.29 


-0.03 


-0.15 




Model 


-0.83 


0.14 


-0.07 




Expt. 


-0.51" 








Expt. 


-0.59' 








Expt. 


-0.68" 








Expt. 


-0.48" 






"Quoted 


in Ref . [TnL 









Hence our prediction for the sign for 3C2 + C3 should be 
sound. In view of the discrepancy with experiment, it 
would be useful to have careful measurements of the rel- 
ative magnitudes of Xmw and > but we are not aware 
of any other work along these lines besides the pioneering 
investigations of Boyd et al^^ 

A convenient measure of the relative importance of 
the two lattice-anharmonicity mechanisms is the ratio 
C2/C3, also included in Table IHTl We expect reason- 
ably accurate ab-initio results for this quantity, since it 
is independent of Xoo i the "weak link" in the calcula- 
tion. Our values clearly cannot be reconciled with those 
of Flytzanis. In the next section we will discuss what this 
implies for the interpretation of the recent experiment of 
Dekorsy et al.^ which attempted to obtain values for the 
parameters C2 and C3 separately for the first time. 



C. Zero-crossings of Xshg the terahertz range 

The quantity |xshg('^)I displayed in Fig.nfor GaAs. 
The dashed line in the upper panel correspond to a sensi- 
ble choice of parameters assembled from the experimen- 
tal and theoretical investigations from the 1960's and 
1970's (these will be referred to as "old parameters"). In 
between the expected strong resonant enhancements at 
a;To/2 and wto, there are two dips, at 5.1 and 7.4 THz, 
the first more pronounced than the second. They result 
from sign reversals of RcxshgC'^) ^'i- O- SHG is 



TABLE 111: Third and fourth columns: parameters C2 and Cs 
combined in a way that relates more directly to experiments. 
Fifth and sixth columns: electro-optic and microwave non- 
linear suceptibilities of Eqs. 112^ and 1)13^ (in pm/V). Rows 
labelled "Model" pertain to the empirical-model calculations 
of Flytzanis. 







3C2 + C3 


C2/C3 


(2) 


(2) 
Xmw 


AlP 


LDA 


-1.68 


-0.026 


28 


-82 




GGA 


-1.66 


-0.022 


26 


-74 


AlAs 


LDA 


-0.83 


-0.028 


50 


-75 




GGA 


-0.77 


-0.023 


48 


-58 


AlSb 


LDA 


-0.23 


0.012 


154 


6 




GGA 


-0.19 


0.017 


145 


25 




Model 


0.93 


-3.33 






GaP 


LDA 


-0.71 


0.146 


88 


-90 




GGA 


-0.78 


0.128 


82 


-69 




Model 


0.27 


-2.22 








Expt. 


0.28" 




20" 


-24" 


GaAs 


LDA 


-0.19 


0.203 


309 


-107 




GGA 


-0.22 


0.175 


240 


-30 




Model 


0.35 


-1.96 








Expt. 


0.39" 




43" 


-51" 



■Refs.Hgl 



observed over a frequency range containing the first zero 
crossing, its frequency loq can be detected as a sharp dip 
in the second harmonic power. If, furthermore, 3C2 + C3 
is known from separate measurements of Xc»\ Xoo^ and 
Xmr^ the remaining free parameter in Eq. (Q), C2/C3, 
can then be adjusted to fit the zero-crossing frequency. 
This was proposed in Ref. as a way of determining C2 
and C3 separately. 

Dekorsy et al^ recently used a free-electron laser to 
measure the far-IR dispersion of the SHG power in GaAs 
from 4.4 to 5.6 THz. They observed the expected reso- 
nance close to WTO /2, followed by a strong drop. Because 
of insufficient filtering of the fundamental signal in the 
detector above 5.6 THz, which masked the SHG signal, 
they were unable to locate precisely the zero crossing, and 
only a lower bound of 5.3 THz was established. Since this 
is slightly above the 5.1 THz predicted for ujq from the 
old parameters, they then discussed how the values of 
C2 and C3 had to be revised to increase ujq to 5.3 THz 
(the assumption being that the lower bound is reason- 
ably close to the actual zero-crossing). They opted to 
leave 3C2 + C3 unchanged at 0.35 (the theoretical value 
from Ref. Hol which is fairly close to the experimental 
0.39); a good fit, shown as a solid line in the upper panel 
of Fig. n was then obtained by changing C2/C3 from 
—2.0 to about —1.23. This amounts to essentially dou- 
bling C3, from —0.07 to —0.14, while changing C2 only 
slightly, from 0.14 to 0.16. 

The dashed line in the lower panel of Fig. ^ shows the 

dispersion obtained with our ab-initio parameters. The 

zero-crossing frequency is raised significantly, to 6.2 THz. 

(2) 

In order to assess the impact of the uncertainty in Xoo 
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GaAs 




Fundamental Frequency (THz) 

FIG. 1: Calculated IR dispersion of |Xshg('^)I 
GaAs, for different choices of the parameters in Eq. Q. 
(Ci,C2,C3) = (-0.59,0.14,-0.07), (-0.59,0.16,-0.13), 
(-0.35,-0.024,-0.12) and (-0.59,-0.041,-0.20) for the 
curves labeled 'Old parameters,' 'Dekorsky,' 'Ab-initio,' and 
'Rescaled ab-initio,' respectively. Following Ref. |^ we set the 
damping parameter 7 to 0.29 THz. The meaning of these 
parameter sets is explained in the text. 



on the dispersion, we show as a solid hne the curve that 
results from reducing Xoo from 472 pm/V to 277 pm/V. 
This affects the C^'s according to Eqs. (jSHSJ, and we 
have chosen the amount of rescaling so as to bring our 
value for Ci into agreement with the experimental num- 
ber from Ref. 34, Ci = —0.59. The zero-crossing fre- 
quency also changes, from 6.2 THz to 5.67 THz, only 
slightly above the measured lower bound of 5.3 THz. 
Clearly, different sets of values for C2 and C3 can lead 
to dispersions with very similar zero crossings (the solid 
lines in the two panels of Fig. , and thus both consis- 
tent with the experimental data, underscoring the need 
for reliable theoretical input. Interestingly, we find that 
for the other HTV compounds no zero crossing occurs for 
tJTo/2 < uj < ojto- This is illustrated in Fig.|21for GaP. 



V. SUMMARY 

We have carried out a detailed ab-initio investigation of 
the IR dispersion of the non-linear susceptibility x*^^^ in 
III-V zincblende semiconductors. The results were com- 
pared with Flytzanis's empirical modeUfi and with ex- 
periment, with particular emphasis on the recent second- 
harmonic generation measurements carried out by Deko- 



GaP 
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Fundamental Frequency (THz) 



FIG. 2: Calculated IR dispersion of IxshgI GaP, using 
the ab-initio parameters without damping (7 = THz). In 
contrast to GaAs (Fig. Ql, no sharp dip is observed between 
the two maxima at u)to/'2 and lotq. 



rsy et alA These authors based the interpretation of their 

data on the parameters obtained in Ref. from model 

calculations. By revising them somewhat, they were able 

(2) 

to obtain a reasonable fit to the IR dispersion of Xshg- 
Instead, we find a completely different set of parameters, 
which however is still consistent with experiment. 

Our findings can be summarized as follows: (i) We 
provide theoretical support to the main qualitative con- 
clusion of Ref.0, that the ratio IC2/C3I between the con- 
tribution from second-order lattice dipole moment (C2) 
versus phonon interaction through the third-order lattice 
potential anharmonicity (C3) is smaller than previously 
thought, (ii) However, we find that this is a consequence 
of not only an increase in | C3 1 )^ but also a significant de- 
crease in IC2I, with the result that the former dominates 
the latter (IC2/C3I -C 1). (iii) The sign of C2 is not 
constant thoroughout the III-V series, and for the two 
most-studied compounds (GaAs and GaP) it is negative, 
contrary to prior understanding, (iv) For all compounds 

except AlSb, we find that the sign of the microwave non- 
(2) 

linear susceptibility Xmw is opposite to that of the optical 

(2) (2) 
iXoo ) and electro-optical (xeo ) coefficients, in agreement 

with early experiments on GaAs and GaPiSiS However, 

our calculated negative sign for 3C2 -I- C3, which is sen- 

(2) (2) 

sitive to the relative magnitudes of Xeo and Xmw, dis- 
agrees with those experiments, (v) The parameter fit 

(2) 

to the IR dispersion of Xshg ^'^li^s on the occurrence 
of a sharp minimum in SHG power between wto /2 and 
loto- The observed zero-crossing in GaAs is reproduced 
by our calculations, but is also consistent with an alter- 
native set of parameters characterized by 3C2 -f C3 > 
and IC2/C3I > 1. For the other compounds considered, 
we find no zero-crossing. 
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APPENDIX A: A METHOD FOR THE DIRECT 
EVALUATION OF Ci AND 3C2 + C3 

As discussed in the main text, by using Eqs. H12() and 
(|13|l one can obtain experimental values for Ci and 3C2 + 

C3 from separate measurements of Xc»\ Xco\ and Xmw- 
In this Appendix we describe the computer-experiment 
analogs of such measurements. 

We define "high-frequency" (optical) fields £ and "low- 
frequency" (static) fields E as follows:-^' when an optical 
field £ is applied the ions are not allowed to move, so that 
the polarization response is purely electronic; in contrast, 
both ions and electrons respond to a static field E, result- 
ing in a total change in polarization that contains ionic 
as well as electronic contributions (but the cell is kept 
strain-free, since we are interested in frequencies above 
the elastic resonances of the medium). In Ref. 13 the 
high-frequency and static dielectric constants, £00 and 
eq, were evaluated by finite differences as first derivatives 
with respect to £ and E, respectively. Extending this to 
second order yields the high-frequency xixl and the low 

,(2) 



frequency (microwave) Xn 

(2) 

Evaluating the electro-optic Xoo requires combining a 
static and an optical field: 



r^(2)i = 

LAeo J 2 dE^ \dSyJ 2 dE^ ^^'^ ^ ^^y' 



(Al) 



Here the action of the static field is described by a total 
derivative as a reminder that the polarization depends on 
a static field both explicitly and implicitly, through the 
atomic positions. Clearly the order of the (partial and 
total) derivatives matters. We can view their combined 
action as a conventional mixed derivative on an auxiliary 
function P(E,£) defined as follows: (i) apply a field E 
and let both electrons and ions respond; (ii) add a field 
£, let the electrons readjust under the total field E + 5 
while keeping the ions fixed in the positions obtained in 
the first step. P(E, £) is defined as the polarization after 
step (ii). Then 



d 
dK 



\^£^ 



dE,d£y 



dSydE, 



(A2) 



E=£=0 ^J^z^'-y 

As before, we apply small fields along [111]: E — e(x 
y-hz), and£ =e(x-hy-hz). Then, defining AP(E, £) 

p(e, e)(x + y + z), we find 

1 d'^p 



X, 



(2) 



which we evaluate as 

d'^p p{e, e) - p{e, ~e) - p{-e, e) + p(-e, -e) 



dede 



4e£ 



+ 0(e^e'). 



(A4) 



We found that a more stringent force tolerance for the 

(2) 

atomic relaxations must be used when evaluating Xmw 
(2) 

than when evaluating Xeo ■ This results from the fact 
that in the latter case the displacements are second- 
order in the field, whereas in the former they are first- 
order Well-converged values were obtained by using a 
force tolerance of 10~^ Ha/bohr for Xmw, while for Xco'' 
10^^ Ha/bohr was sufficient. 

The values for Xeo'' and Xmw obtained using this 
method agree to within 1 pm/V with the ones obtained 
from the separate calculation of Ci , C2 and C3 described 
in Sec, mil This provides an internal consistency check of 
our calculations. Although it is somewhat more expen- 
sive and does not provide as much information (e.g., it 
does not produce separate values for C2 and C3), the 
method described in this Appendix provides a simple 
means of controlling the mechanical boundary conditions 
under applied fields. Although we only considered atomic 
displacements, the same strategy can be extended to 
strain deformations. In this way one can easily compute, 
for example, the clamped (strain-free) and undamped 
(stress-free) electro-optic coefficients, or the static x'^' 
including the strain response^ 



APPENDIX B: IMPROVING THE 
CONVERGENCE WITH RESPECT TO fc-POINT 
SAMPLING 

Although total-energy ground state calculations for 
insulators converge exponentially fast with respect to 
fc-point sampling, for finite-field calculations the con- 
vergence is considerably slower This results from 
the discretized Berry-phase (DBF) polarization expres- 
sion (Eq. (|B2|I below) used in the electric enthalpy 
functional and as a consequence the second field- 
derivatives in Eqs. l|HJl and (|5Jl also converge slowly. In 
order to alleviate this problem we used in our finite-field 
calculations a modified DBF expression for the polariza- 
tion, Eq. HB4|I below. 

For notational simplicity we limit our discussion to the 
case of a single valence band in one dimension. The 
electronic polarization of a bulk insulator under periodic 
boundary conditions can be written, by analogy with the 
dipole moment of a molecule, as Pei = —e{x)/L, where 
L is the length of the periodic box. The Berry-phase 
expression for (x) ia^ 



2-K 



Im ln(vl'|e'^|\l'). 



(Bl) 



A dede' 



(A3) 



where is a many-body insulating wavefunction. For 
a finite system in a supercell, evaluating the dipole 
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moment with this expression amounts to replacing the 
non-periodic operator x with the periodic operator 
(L/27r) sin(27ra;/Z/)4ii The difference between the two 
near the origin, where the molecule is located, is of order 
Using for ^ a Slater determinant of single-particle 
Bloch states, we recover from this expression the King- 
Smith- Vanderbilt DBP expression for the polarization of 
a band insulator ;.2£ 



Pel = Im In TT det S{k,, fc,+i) + ©(l/L^), (B2) 



N-1 



where S{ks,ks+i) = (ufcjufe,+i)- 
sion for (x) is 



of predicting the converged values, with errors of usually 
around 1%. This procedure requires several calculations 
at different N, starting at = 12 to avoid the contri- 
bution from higher-order terms in The modified 
polarization expression produces results of similar accu- 
racy with a calculation for a single value of N. We found 
that iV = 16 usually provides accurate values for the cal- 
culation of nonlinear susceptibilities to within 1 pm/V. 
As we can see from the graph, //^^^ calculated in this 
way converges rapidly: The resulting value at A'^ = 6 is 
closer to the converged value of 0.284 nC/m than that 
calculated with the conventional functional at = 20. 

(2) 

An alternative expres- ^^'^ improvement for aro and Xoo is less dramatic, but 



ZTT 



-Im ln(\I>|e'" 



I*) - -Im ln(*le*^|«') 



(B3) 

which is correct to ©(l/L**), as can be seen using the 
same type of reasoning as in Ref. ^3 Eq. (|B3p leads to 
a modified DBP polarization formula, 



e 

2^ 



-Im In Y[ det S{ks,ks+i) 



s=0 



-Im In Yl det S{ks,k, 



s+2) 



s=0 



0{1/L^). 
(B4) 



This is closely related to the expression obtained by 
combining Richardson's extrapolation, Eq. (|20|l . with 
Eq. (|B2|I . viewed as a finite-difference representation of 
the fc-derivative in the continuum Berry-phase formula.**^ 
We show in Figure |21 the convergence with the size 
N X N X N of the shifted Monkhorst-Pack grid of the 
second-order quantities Xm"* , axo , and /x'^^^ (these are the 
third derivatives of the electric enthalpy which involve at 
least one field derivative) . For each quantity, we plot as a 
function of 1/A^^ the value obtained using both the con- 
ventional polarization expression HB2() (dashed line) and 
the modified expression IIB4|) (solid line) . Notice that the 

values for ■: c^tO: a-nd m'^'' obtained from the former 
fall nearly on a straight line. Extrapolating to = cxo 
through a least squares fit against 1/A^^ is a reliable way 
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FIG. 3: Comparison of the convergence of , Qto, and 
/x'^-* for AlP with respect to the fc-point mesh when using the 
conventional discretized Berry-phase expression <B2ll (dashed 
lines) versus using the modified expression <B4ll (solid lines). 



in general there is a clear improvement for most materi- 
als. The notable exception is GaAs, the material with the 
smallest gap, where there is an improvement for all quan- 

(2) 

titles except Xoc i which actually converges more slowly 
when using the modified fourth-order expression. We 
have found empirically that this expression works best 
for large band gap materials. 
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